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a first signal representing physical characteristics 
of a reservoir, obtaining a residual vector, (Ro) 
from the first signal (representing errors associated 
with a system of nonlinear equations describing 
the reservoir), and a first matrix, (Ao), (78), 
(representing the sensitivity of the residual vector 
to changes in a system of nonlinear equations), 
recursively decomposing matrix (Ao), into a lower 
block triangular matrix, an upper block triangular 
matrix, and a diagonal matrix, and generating, 
(80), a second matrix (Mo) that is an approximation 
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AN IMPROVED SIMULATION METHOD AND APPARATUS 

5 BACKGROUND OF THE INVENTION 

The subject matter of the present invention relates to a simulator software method 
and apparatus, embodied in an earth formation reservoir simulator, for simulating an 
earth formation reservoir containing liquids and/or gases by solving a system of 
10 linear equations that characterize physical aspects of an oil and/or gas field, the 
amount of time required by the reservoir simulator to solve the system of linear 
equations and to thereby simulate the earth formation reservoir being reduced 
relative to prior art simulation methods practiced by prior art simulators. 

15 Oil and gas is produced from underground rock formations. These rocks are porous, 
just like a sponge, and they are filled with fluid, usually water. This porous 
characteristic of rocks is known as porosity. These rocks in addition to being porous 
have the ability to allow fluids to flow through the pores, a characteristic measured 
by a property called permeability. When oil (or gas) is trapped in such formations, it 

20 may be possible to extract it by drilling wells that tap into the formation. As long as 
the pressure in the well is lower than that in the rock formation, the fluids contained 
in the pores will flow into the well. These fluids may then flow naturally up the well 
to the surface, or the flow up the well may have to be assisted by pumps. The 
relative amounts of oil, gas and water that are produced at the surface will depend on 

25 the fraction of the rock pore space that is occupied by each type of fluid. Water is 
always present in the poies, but it will not flow unless its volume fraction exceeds a 
threshold value that varies from one type of rock to another. Similarly, oil and gas 
will only flow as long as their volume fractions exceed their own thresholds. 

30 The characteristics of the rock (including porosity and permeability) in an oil 

reservoir vary greatly from one location to another. As a result, the relative amounts 
of oil, gas and water that can be produced will also vary from reservoir to reservoir. 
These variations make it difficult to simply predict the amount of fluids and gases a 
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reservoir will produce and the amount of resources it will require to produce from a 
. particular reservoir. However, the parties interested in producing from a reservoir 
need to project the production of the reservoir with some accuracy in order to 
determine the feasibility of producing from that reservoir. Therefore, in order to 
5 accurately forecast production rates from all of the wells in a reservoir, it is 
necessary to build a detailed mathematical model of the reservoir's geology and 
geometry. 

A large amount of research has been focused on the development of reservoir 
10 simulation tools. These tools include mathematical and computer models that 

describe and which are used to predict, the multiphase flow of oil and gas within a 
three dimensional underground formation (a'field*). Reservoir tools use empirically 
acquired data to describe a field. These data are combined with and manipulated by 
mathematical models whose output describes specified characteristics of the field at 
15 a future time and in terms of measurable quantities such as the production or 
injection rates of individual wells and groups of wells, the bottom hole or tubing 
head pressure at each well, and the distribution of pressure and fluid phases within 
the reservoir. 



20 



The mathematical model of a reservoir is typically done by dividing the reservoir 
volume into a large number of interconnected cells and estimating the average 
permeability, porosity and other rock properties for each cell. This process makes 
use of seismic data, well logs, and rock cores recovered when wells are drilled. 
Production from the reservoir can then be mathematically modeled by numerically 
25 solving a system of three or more nonlinear, partial differential equations describing 
fluid flow in the reservoir. 

Computer analysis of production from an oil reservoir is usually divided into two 
phases, history matching and prediction. In the history matching phase, the past 
30 production behavior of the reservoir and its wells is repeatedly modeled, beginning 
with initial production and continuing up to the present time. The first computer run 
is based on a geological model as described above. After each run, the computer 
results are compared in detail with data gathered in the oil field during the entire 
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period of production. Geoscientists modify the geological model of the reservoir on 
the basis of the differences between computed and actual production performance 
and rerun the computer model. This process continues until the mathematical 
reservoir model behaves like the real oil reservoir. 



Once a suitable history match has been obtained, production from the oil reservoir 
can be predicted far into the future (sometimes for as long as 50 years). Oil recovery 
can be maximized and production costs minimized by comparing many alternative 
operating plans, each requiring a new run of the computer model. After a field 
10 development plan is put into action, the reservoir model may be periodically rerun 
and further tuned to improve its ability to match newly gathered production data. 

When sufficient data is obtained about the reservoir, characteristics of a reservoir 
can be mathematically modeled to predict production rates from wells in that 
15 reservoir. The gross characteristics of the field include the porosity and 

permeability of the reservoir rocks, the thickness of the geological zones, the 
location and characteristics of geological faults, relative permeability and capillary 
pressure functions and such characteristics of the reservoir fluids as density, 
viscosity and phase equilibrum relationships. From this data, a set of continuous 
20 partial differential equations (PDEs) are generated that describe the behavior of the 
field as a function of time and production parameters. These production parameters 
include the locations of wells, the characteristics of the weirs completions, and the 
operating constraints applied to the wells. Operating constraints may include such 
as the production rate of a particular fluid phase, the bottom hole pressure, the 
25 tubing head pressure, or the combined flow rates of a group of wells. These 

constraints may be applied directly by data or by means of another simulator that 
models the flow of fluids in the surface equipment used to transport the fluids 
produced from or injected into the wells. However, because only the simplest 
system of PDEs can be solved using classic or closed-form techniques (e.g., a 
30 homogeneous field having circular boundaries), a models PDEs are converted into a 
set of non-linear approximations which are then solved numerically. One 
approximation technique is the finite difference method. In the finite difference 
method, reservoir PDEs are converted into a series of difference quotients which 
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divide a reservoir into a collection of discrete three dimensional cells, which are then 
solved for at discrete times to determine (or predict) the value of reservoir 
characteristics such as pressure, permeability, fluid fractions, and at a later time. 

5 Within the computerized reservoir simulator, reservoir performance is modeled in 
discrete increments of time. Each so-called timestep advances the solution from a 
previous point in time, where all variables are known, to a future point in time, 
where all variables are unknown. This process is repeated until the entire time 
period of interest has been modeled. Within each timestep it is necessary to solve a 

10 huge system of nonlinear equations that models fluid flow from cell to cell and 

through the wells. (With current technology it is possible to include several million 
cells in the reservoir model.) Solutions to the system of nonlinear equations are 
obtained by Newton iteration. In each such iteration the system of nonlinear 
equations is approximated by a system of linear equations, which must be solved by 

15 yet another iterative procedure. 

A general outline of the operation of a reservoir simulator follows (refer to 
figure 8a). Reservoir data and rock core data are used to describe a computational 
grid and the properties of the reservoir rocks. This is combined with data on the 
20 physical properties of the fluids contained in the reservoir and used to compute the 
initial distributions of pressure and fluid saturations (volume fractions) as well as the 
composition of each fluid phase. Time varying data, such as the locations and 
characteristics of wells, production and injection flow rate controls, and simulator 
control information, is read from a data base. Using the current pressure, saturation, 
25 and fluid compositions for each grid cell, the partial differential equations describing 
mass balances are approximated by finite differences resulting in two or more 
nonlinear algebraic equations for each cell. In addition, these nonlinear equations 
are linearized by means of Newton's method. The resulting system of linear 
equations is solved iteratively using methods described in this specification. After 
30 the linear equations have been solved, a test is performed to determine whether all of 
the nonlinear terms in the finite difference equations have converged. If not, the 
simulator returns to a previous step wherein the partial differential equations 
describing mass balances are approximated by finite differences. However, if the 
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nonlinear terms have converged, the simulator updates values to complete the 
current timestep. Then, the simulator tests to determine whether the desired ending 
time in the simulation has been reached. If not, the simulator returns to a previous 
step to read new time varying data and begins the next timestep. However, if the 
5 endpoint of the simulation has been reached, then, the simulator completes output 
operations and the run is finished. 

As reservoir simulations grow in complexity (e.g., number of parameters) and size 
(e.g., number of cells), solving the resulting system of linear equations, represented 
10 by the matrix equation A „ requires an increasingly large effort. For example, 
the work involving the iterative solution increases with the square of the number of 
parameters per cell. In many reservoir models, the time required to solve systems of 
linear equations (step 56 in figure 8a) can be a limiting factor on the simulation's 
usefulness. 

15 

In connection with reservoir modeling, well logging operations are performed in the 
formation thereby producing well log data, and seismic operations are performed on 
the formation thereby producing seismic data. The seismic data is reduced thereby 
producing reduced seismic data. The well log data and the reduced seismic data are 
20 introduced, as input data, to a computer workstation which stores a gridding 

software and a simulator software. A gridding software, hereinafter known as'the 
Flogrid software" or the'Flogrid gridding software", is disclosed in prior pending U.S. 
patent application serial number 09/034,701, filed in the U.S. on March 4, 1998, 
which is based on a Great Britain patent application number 9727288.4 filed 
25 December 24, 1997, the disclosure of which is incorporated by reference into this 
specification. The"Flogrid' gridding software includes another gridding software 
known as'Tetragrid'. The'Petragrid' gridding software is disclosed in prior pending 
U.S. patent application serial number 08/873,234 filed June 1 1, 1997, the disclosure 
of which is also incorporated by reference into this specification. The gridding 
30 software will respond to the reduced seismic data and the well log data by gridding 
the earth formation which was subjected to the well log operation and the seismic 
operation. The type of grids imposed on the earth formation include structured 
(approximately rectangular) grids and unstructured (tetrahedral) grids. A property, 
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such as permeability or water saturation, is assigned to each cell or grid block of the 
grid. As a result, a set of output data is generated by the gridding software, the set 
of output data including the plurality of cells/grid blocks of the grid and the 
respective plurality of properties associated with each of the cells of the grid. 

5 

The set of output data from the gridding software are introduced, as input data, to a 
reservoir simulator software. The reservoir simulator software will respond to the 
set of output data from the gridding software by generating a plurality of simulation 
results which are associated, respectively, with the plurality of cells/grid blocks of 
10 the grid received from the gridding software. The plurality of simulation results 

and the plurality of cells/grid blocks associated therewith, generated by the reservoir 
simulator software, will be displayed on a 3D viewer of the workstation for 
observation by a workstation operator. Alternatively, the plurality of simulation 
results and the plurality of cells/grid blocks associated therewith can be recorded for 
15 observation by a workstation recorder. 

The reservoir simulator software can model an oilfield reservoir. For example, in 
the Society of Petroleum Engineers (SPE) publication number 28545, concerning a 
transient tool for multiphase pipeline and well simulation, dated 1994, the authors 
20 have solved for pressure losses along a single pipeline using a technique related to 
conservation of material and conservation of pressure. A similar technique has 
been applied to a network of pipelines or flowlines in the Society of Petroleum 
Engineers (SPE) publication number 29125, authored by Litvak and Darlow. In this 
publication, the authors (Litvak and Darlow) have taken a network model (i.e., a 

25 network of pipelines) in which the pressure losses along the network branches can 
either be calculated from tables or from an analytical model, and the analytical 
model solves for three (3) conservations and pressures. In addition, in an article by 
the'Society of Petroleum Engineers" (SPE) 12259, each well being modeled in that 
article was characterized by three (3) variables: pressure, water fraction, and gas 

30 fraction. 

As noted earlier, in many reservoir models, the time required to solve systems of 
linear equations (step 56 in figure 8a) can be a limiting factor on the simulation's 
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usefulness. Accordingly, it would be beneficial to have a method that can be used to 
efficiently solve large systems of linear equations. It is the solution of this system of 
linear equations that is the subject of the invention of this specification. 

5 Accordingly, a new and improved reservoir simulator software is needed wherein 
the amount of time required to execute the new and improved reservoir simulator 
software is much less than the amount of time required to execute the reservoir 
simulator software of the prior art. 

10 SUMMARY OF THE INVENTION 

Accordingly, it is a primary object of the present invention to provide an improved 
method of performing computer simulations. 

15 It is a primary aspect of this invention to provide a simulation program that includes 
a'linear Solvef 'method and apparatus which operates with reduced processing time 
relative to prior art simulators. 

It is a further aspect of this invention to provide a simulation program that includes a 
20 'linear Solvef 'method and apparatus which practices a method for solving systems 
of differential equations having substantially fewer steps than conventional methods 
of solving differential equations. 

It is a further aspect of this invention to provide a simulation program which 
25 includes a "Linear Solver" method and apparatus that generates a set of simulation 
results prior to completely solving a system of differential equations. 

In accordance with the primary object and other aspects of the present invention, an 
improved "Linear Solver" method and apparatus is disclosed for simulating 
30 conditions and characteristics that are mathematically modeled using linear and non- 
linear systems of differential equations. These systems of differential equations are 
very large and complex and require substantial processing time to solve. The 
processing time required to execute a simulation program can have an effect on the 



7 



WO 00/75854 



PCT/USOO/15056 



accuracy of the simulation results. The present invention provides a faster method 
to solve differential equation systems by generating a smaller approximated matrix 
that is an approximation of the system of equations. This smaller approximated 
matrix is smaller than another matrix that represents the original system of 

5 equations. From the smaller approximated matrix, information that is usually 
determined from the solved equations is obtained prior to the step of solving the 
differential equation system. The reduced size of the smaller approximated matrix, 
coupled with the fact that some of the final information is obtained prior to 
completely solving the system of differential equations, results in a faster simulation 

10 method. 

Consequently, the'linear Solvef 'of the present invention provides a method for 
solving a system of linear equations representing physical characteristics of an oil 
and/or gas reservoir. This method includes receiving a first signal representing 
15 physical characteristics of the reservoir, obtaining a residual vector f o (representing 
errors associated with a system of nonlinear equations describing the reservoir) and 
a first matrix Ao (representing the sensitivity of the residual vector f o to changes 
in the system of nonlinear equations) from the received first signal; recursively 
decomposing matrix Ao in» a lower block triangular matrix, an upper block 
20 triangular matrix, and a diagonal matrix; and generating a second matrix Mo that 
is :an approximation to matrix Ao- The invention further provides generating an 
approximate solution for the system of linear equations (in the form of a signal) 
without directly computing the matrix multiplication A 0 * = * required of 
conventional techniques. 



25 



At each time step of the simulation, it is necessary to solve a system of linear 
equations represented by A** - TV A> « a matrix obtained from the data input 
into the simulation program. Ao is a SP 8 ** matrix in that there are few non-zero 
elements in the matrix. Because of the size of the matrix, it would be very costly to 
30 solve the system of differential equations associated with this matrix Ao directly. 
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Therefore, Ao is decomposed into smaller matrix components that will be easier to 
factor, including an upper matrix component, a lower matrix component and 
diagonal matrix components. 

5 An approximate matrix A/ 0 is generated from the matrix components. This matrix 
M 0 can be factored. Once the factorization process is complete, this matrix M o 
can be solved to get a result similar to the result obtained from solving the original 
matrix A 0 . However, solving the M o matrix requires much less time than the 
original matrix^ would require. This process is possible because data normally 

10 determined from solving the original matrix Ao is obtained from information 
determined during the factorization of matrix Mo • Therefore, the process time 
required to solve the matrix Mo is much less because some of the information 
obtained from solving the original matrix Ao has already been determined during 
the factorization of matrix Mo • 

15 

The inventive concept may be embodied in a computer program stored in any media 
that is readable and executable by a computer system. Illustrative media include 
magnetic and optical storage devices. Illustrative programming languages include 
Fortran and C Illustrative computer systems include personal computers, 
20 engineering workstations, minicomputers, mainframe computers, and specially 
designed state machines. 

Further scope of applicability of the present invention will become apparent from 
the detailed description presented hereinafter. It should be understood, however, 
25 that the detailed description and the specific examples, while representing a 

preferred embodiment of the present invention, are given by way of illustration only, 
since various changes and modifications within the spirit and scope of the invention 
will become obvious to one skilled in the an from a reading of the following 
detailed description. 
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BRTPF DESCRIPTION OF TH F DRAWINGS 

A full understanding of the present invention will be obtained from the detailed 
description of the preferred embodiment presented hereinbelow, and the 
5 accompanying drawings, which are given by way of illustration only and are not 
intended to be limitative of the present invention, and wherein: 

figure 1 illustrates a seismic operation for producing a reduced seismic data output 
record, the seismic operation of figure 1 including a data reduction operation; 

10 

figure 2 illustrates a wellbore operation for producing a well log output record; 

figure 3 illustrates a computer system for performing the data reduction operation of 
figure 1; 

15 

figures 4 and 5 illustrate a workstation adapted for storing a"Flogrid"software and an 
•Eclipse" simulator software; 

figures 6 and 6a illustrate a more detailed construction of the"Flogrid"software of 
figure 5 which is adapted for generating output data for use by the'Eclipse" simulator 
software, the Eclipse simulator software including a'linear Solver method and 
apparatus" in accordance with the present in vention; 

figure 7 illustrates an example of a typical output display generated by the'Eclipse" 
simulator software of figure 6 which is displayed on the 3D viewer of 
figure 6; 

figure 8a illustrates a prior art approach or method of performing reservoir 
simulation which has been practiced by prior art reservoir simulators; 

figure 8b illustrates the Eclipse simulator software of figure 5 including the Linear 
Solver of the present invention; 
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figures 9 through 13 illustrate the 4 l-inear Solve*' method and apparatus in accordance 
with the present invention which is embodied in the'Eclipse simulator software of 
figures 5 and 6, and wherein: 

figure 9 illustrates the approach or method of performing reservoir simulation 
practiced by the'Eclipsd' simulator software of present invention; 

figure 10 outlines a technique or method for solving a system of linear equations; 

figure 11 illustrates a matrix decomposition technique or method practiced by the 
'Eclipsd* simulator software of the present invention; 

figure 12 illustrates a factorization technique practiced by the'Eclipsd' simulator 
software in the present invention; and 

figure 13 illustrates a method for producing an approximate solution to an original 
system of linear equations in accordance with an inventive concept of the present 
invention. 

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT 

The'Eclipsd' simulator software of the present invention receives output data from 
the'Flogrid' simulation gridding software and, responsive thereto, the'Eclipsd* 
simulator software generates a set of simulation results which are displayed on a 3D 
viewer. The'Eclipsd' simulator software includes the'linear Solvef method and 
apparatus in accordance with the present invention. The'linear Solved apparatus of 
the present invention includes a special apparatus which allows the simulation 
processing time of the'Eclipsg' simulator to be substantially less than the processing 
time associated with the prior art reservoir simulators. 

This specification is divided into two parts: (1) a Background discussion which 
provides background information relating to the performance of a seismic operation 
and a well logging operation adapted for generating seismic and well logging data, 
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the seismic and well logging data being provided as input data to a workstation that 
stores a"Flogrid'simulation gridding software and an"Eclipse" simulator software, and 
(2) a more detailed description of the'Eclipse" simulator software, the"Eclipse" 
simulator software including a'linear Solvef'method and apparatus in accordance 
with the present invention which causes or allows the simulation processing time of 
the Eclipse simulator to be substantially less than the simulation processing time 
associated with other prior art simulators. 

Background Discussion 
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Referring to figure 1, a method and apparatus for performing a seismic operation is 
illustrated. During a seismic operation, a source of acoustic energy or sound 
vibrations 10, such as an explosive energy source 10, produces a plurality of sound 
vibrations. In figure 1 , one such sound vibration 12 reflects off a plurality of 
15 horizons 14 in an earth formation 1 6. The sound vibration(s) 12 is (are) received in 
a plurality of geophone-receivers 18 situated on the earth's surface, and the 
geophones 18 produce electrical output signals, referred to as"data received'20 in 
figure 1, in response to the received sound vibration(s) 12 representative of different 
parameters (such as amplitude and/or frequency) of the sound vibration(s) 12. The 
20 'data received*20 is provided as"input data'to a computer 22a of a recording truck 22, 
and, responsive to the'mput data", the recording truck computer 22a generates a 
'teismic data output record'24. Later in the processing of the seismic data output 
record 24, such seismic data undergoes "data reductiori'30 in a mainframe computer, 
and a'teduced seismic data output record'24a is generated from that data reduction 
25 operation 30. 

Referring to figure 2, a well logging operation is illustrated. During the well logging 
operation, a well logging tool 34 is lowered into the earth formation 1 6 of figure 1 
which is penetrated by a borehole 36. In response to the well logging operation, 
30 well log data 38 is generated from the well logging tool 34, the well log data 38 
being provided as'lnput data" to a computer 40a of a well logging truck 40. 
Responsive to the well log data 38, the well logging truck computer 40a produces a 
'well log output record' 42. 
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Referring to figure 3, the seismic data output record 24 of figure 1 is provided as 
'input data" to a mainframe computer 30 where the data reduction operation 30 of 
figure 1 is performed. A mainframe processor 30a will execute a data reduction 
software 30b stored in a mainframe storage 30b. When the execution of the data 
reduction software 30b is complete, the reduced seismic data output record 24a of 
figures 1 and 3 is generated. 

Referring to figures 4 and 5, a workstation 44 is illustrated in figure 4. A storage 
0 medium 46, such as a CD-Rom 46, stores .software, and that software can be loaded 
into the workstation 44 for storage in the memory of the workstation. In figure 5, 
the workstation 44 includes a workstation memory 44a, the software stored on the 
storage medium (CD-Rom) 46 being loaded into the workstation 44 and stored in 
the workstation memory 44a. A workstation processor 44d will execute the 
15 software stored in the workstation memory 44a in response to certain input data 
provided to the workstation processor 44d, and then the processor 44d will display 
or record the results of that processing on the workstation'tecorder or display or 3D 
viewef'44e. The input data, that is provided to the workstation 44 in figure 5, 
includes the well log output record 42 and the reduced seismic data output record 
20 24a. The"well log output record'42 represents the well log data generated during the 
well logging operation in an earth formation of figure 2, and the'Yeduced seismic 
data output record'24a represents data-reduced seismic data generated by the 
mainframe computer 30 in figure 3 in response to the seismic operation illustrated in 
figure 1. In figure 5, the software stored on the storage medium (CD-Rom) 46 in 
25 figure 5 includes a'Tlogrid' software 46a and an"Eclipse" simulator software 46b. 

When the storage medium (CD-Rom) 46 is inserted into the workstation 44 of figure 
5, the'Flogrid' software 46a and the'Eclipse" simulator software 46b, stored on the 
CD-Rom 46, are both loaded into the workstation 44 and stored in the workstation 
memory 44a. The'TFlogrid' software 46a is fully described and set forth in pending 
30 U.S. patent application serial number 09/034,701 , filed in the U.S. on March 4, 
1998, which is based on prior pending Great Britain patent application number 
9727288.4 filed December 24, 1997, the disclosure of which is incorporated by 
reference into this specification. When the workstation processor 44d executes the 
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Flogrid software 46a and the Eclipse simulator software 46b, the'Eclipse" simulator 
software 46b responds to a set of more accurate grid cell property information 
associated with a respective set of grid blocks of a structured simulation grid 
generated by the'Flogrid' software 46a by further generating a set of more accurate 
5 simulation results which are associated, respectively, with the set of grid blocks of 
the simulation grid Those simulation results are displayed on the 3D viewer 44e of 
figure 5 and can be recorded on a recorder 44e. 

Referring to figures 6 and 6a, in figure 6, the Flogrid software 46a and the Eclipse 
10 simulator software 46b are illustrated as being stored in the workstation memory 44a 
of figure 5 In addition, the"simulation results" 48, which are output from the Eclipse 
simulator software 46b, are illustrated as being received by and displayed on the 3D 
viewer 44e. The Flogrid software 46a includes a reservoir data store, a reservoir 
framework, a structured gridder, an unstructured gridder, and an upscaler, all of 
15 which are fully discussed in the above referenced prior pending U.S. patent 
application serial number 09/034,701, filed in the U.S. on March 4, 1998, the 
disclosure of which has already been incorporated by reference into this 
specification. In figure 6, a set of'simulation grids and properties associated with 
the grids"47, generated by the Upscaler and the'Petragrid'unstructured gridder, are 
20 received by the Eclipse simulator software 46b. In response, the Eclipse simulator 
software 46b generates a"set of simulation results associated, respectively, with a set 
of grid blocks of the simulation grids"48, and the simulation results and the 
associated grid blocks 48 are displayed on the 3D viewer 44e. 

25 In figure 6a, the Flogrid software 46a generates a set of output data 47 comprising a 
plurality of grid cells and certain properties associated with those grid cells. That 
output data 47 is provided as input data to the Eclipse simulator software 46b. Some 
other programs 49 provide other input data to the Eclipse simulator software 46b. In 
response to the output data 47 (comprised of a gridded earth formation including a 
30 plurality of grid cells and certain properties associated with each of the grid cells), as 
well as the other output data from the other programs 49, the Eclipse simulator 
software 46b generates a set of 'simulation results*^, the simulation results 48 
including the plurality of grid cells and a plurality of simulation results associated. 
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respectively, with the plurality of grid cells. The aforementioned plurality of grid 
cells and the plurality of simulation results associated, respectively, with the 
plurality of grid cells are displayed on the 3D Viewer 44e of figure 6 and 6a. 

Referring to figure 7, an example of the simulation results 48 (i.e., the "jiurality of 
grid cells and the plurality of simulation results associated, respectively, with the 
plurality of grid cellsT48) which are displayed on the 3D viewer 44e of figures 5 and 
6 and and 6a, is illustrated in figure 7. 

Detailed Description of the Eclipse simulato r software 46b of figures 5 and 6 
including the "Linear Solver" method and apparatus in accordance with the present 
invention 

Referring to figure 8a, a general outline of the operation of a prior art reservoir 
simulator is discussed below with reference to figure 8a. In figure 8a, reservoir data 
42 and 24a of figure 5 and rock core data are used to describe a computational grid 
and the properties of the reservoir rocks. This data is combined with data relating to 
the physical properties of the fluids contained in the reservoir, the combined data 
being used to compute the initial distributions of pressure and fluid saturations 
(volume fractions) as well as the composition of each fluid phase, block 50 in figure 
8a, Time varying data, such as the locations and characteristics of wells, production 
and injection flow rate controls, and simulator control information is read from a 
data base, block 52. Using the current pressure, saturation, and fluid compositions 
for each grid cell, the partial differential equations describing mass balances are 
approximated by finite differences in block 54 which results in two or more 
nonlinear algebraic equations for each grid cell. Also, in block 54, these nonlinear 
equations are linearized by means of Newton's method. In block 56, the resulting 
system of linear equations is solved iteratively, using methods described in this 
specification. After the linear equations have been solved, there is a test in block 58 
to determine whether all of the nonlinear terms in the finite difference equations 
have converged. If not, the simulator returns to block 54. If the nonlinear terms in 
the finite difference equations have converged, the simulator moves to block 60 and 
updates values to complete the current timestep. In block 62, the simulator tests to 
determine whether the desired ending time (i.e., the stop time) in the simulation has 
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been reached. If not, the simulator advances time to the next level, block 64, and 
then it returns to block 52 to read new time varying data and to begin the next 
timestep. If the endpoint of the simulation has been reached, then, the simulator 
completes output operations and the run is finished, block 66. 

Referring to figure 8b, as noted earlier, the Eclipse simulator software 46b of figures 
5 and 6 includes the "Linear Solver" 70 in accordance with the present invention. 
The Linear Solver 70 substantially reduces the simulation processing time which is 
normally associated with prior art reservoir simulators. 

As noted earlier with reference to figure 8a, blocks 54 and 56, nonlinear partial 
differential equations are linearized (block 54 of figure 8a) and the resulting system 
of linear equations is solved iterative!* using methods discussed in this specification 
(block 56 of figure 8a). If the nonlinear terms are not converged (block 58 in figure 
8a), the following discussion with reference to figures 9 through 13 will discuss 
those methods for solving a system of non-linear partial differential equations. 



Referring to figure 9, a description of the Linear Solver 70 of the present invention 
of figure 8b will be discussed below in the context of a novel method for solving a 
20 system of nonlinear partial differential equations that represent an oil or gas 
reservoir, that novel method being illustrated in figure 9. 

The basic distinction between this invention and the conventional simulation 
methods is the way step 56 of figure 8a is performed (i.e.,«lteratively solve the linear 
25 equations'). Systems of differential equations are very large and to solve them 

directly requires solving the equation Ao* = B for each time ste P of the simulation 
process. The conventional method would require substantial time and effort To 
reduce the amount of processing time, the conventional simulation methods would 
require a reduction the amount of input data, which prevents a more accurate 
30 simulation. However, the simulation method (and apparatus) of the present 
invention performed by the Linerar Solver 70 performs steps that enable one to 
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directly avoid the equation A 0 x = B by solving an approximate system of equations 
in less time and without sacrificing accuracy in the results. 

In figure 9, the Linear Solver 70 simulation method of present invention: 
(1) converts the system of non-linear equations into a system of linear equations, 
block 76 of figure 9; (2) solves the system of linear equations; and then (3) solves 
the system of non-linear equations using the solved system of linear equations 
obtained during step (2). 

The method practiced by the Linear Solver 70 of figure 8b of the present invention, 
for solving the system of linear equations from step (2) above, involves four basic 
steps that distinguishes the simulation method of the present invention from the 
conventional simulation methods. These four basic steps are set forth below with 
reference to figure 9 as follows: 

(1) Obtain a matrix Ao of tne system of linear differential equations, block 78 in 
figure 9; 

(2) Decompose this matrix ^into components [i.e., an upper block component, a 
lower block component, and a diagonal component], to facilitate factorization of 
the matrix Ao> until a new factorable matrix Mo is obtained, block 80 in figure 
9; 

(3) Compute a factorization of this new factorable matrix M 0 , block 82 in figure 
9, to produce a system of linear equations that approximates the original matrix 

(4) Solve the system of linear equations from step (3) which are represented by 
matrix Mo> block 84 of figure 9. 

This solution to matrix Mo gi ves an approximate but accurate solution to the 
original matrix Ao- This solution to matrix Mo is ™ en us ed to solve the system of 
non-linear equations, block 86 in figure 9. The advantage of solving the 
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matrix Mo is that this matrix Mo requires much less time to solve. In addition, 
the accuracy of the matrix Mo « not a compromise because much of the 
Mnformatiort'that is usually determined from calculating the matrix A 0 " calculated 
during the factorization step of matrix Mohock 82 of fig. 9). Because this 
•Informatiori'is already calculated from the factorization step of matrix M o • block 
82, there is no need to calculate this"informatiori' when solving matrix M o . Wock 
84. Therefore it takes less time to solve this matrix Mo because there are fewer 
components to calculate. 

In figure 9, generate a system of approximate linear equations from the original 
system of nonlinear equations, block 76 of figure 9, and generate a matrix Ao from 
these linear equations, block 78. The matrix Ao . generated from the linear 
equations in block 78 of figure 9, is sparse in that it contains few non-zero elements. 
The small number of non-zero elements in this matrix A 0 and the fact that this 
matrix Ao ^ very large and complex make it very difficult to solve the resulting 
system of equations directly. In order to solve this matrix A 0 * *» matrix A * 
first decomposed into components (upper block, lower block and a diagonal), block 
80 of figure 9. This decomposition divides the matrix Ao ™ to a P luralitv of smaller 
matrix blocks. Each of the plurality of smaller matrix blocks of the matrix Ao * 
, then recursively decomposed (the upper block from bottom to top and the lower 
block from top to bottom)) until aresulting matrix Mo* generated which can te 
efficiently factored exactly (block 80). This decomposition of the matrix Ao results 
inanew matrix Mo that is an approximation of the original matrix A 0 - The next 
step, block 82 of figure 9, involves the computation of a factorization of the matrix 
:5 Mo^hecomputationofthefactorizationofthematrix Mo generating a system of 
equations which is an approximation of the original set of linear equations and 
which can be more efficiently solved. This factorization step (computing the 
factorization of matrix M . . block 82 of figure 9) generates the following equation: 
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The solution of the above equation"M 0 = (Lo M \ + D( Mi + i/ 0 ) " for 
matrix yV</o» block 84 of fi g ure 9 > can be used to determine the solution to the 
original system of non-linear equations, block 86 of figure 9. 

As indicated in block 82 of figure 9 (and the above step 3), one widely used 
technique to solve systems of linear equations is'Viested factorizatiorP which is 
defined to be an iterative technique for determining an approximate matrix 
corresponding to the system's coefficient matrix. A system of linear equations 
representing reservoir characteristics may be solved by exploiting certain, and 
heretofore unrecognized, characteristics of the nested factorization method. The 
following embodiment of this inventive concept, relating to oil and gas reservoir 
modeling, is illustrative only and is not to be considered limiting in any respect. 

Referring to figures 10 through 13, these figures will discuss, in detail, some of the 
individual steps of figure 9 which are practiced by the Linear Solver 70 of the 
present invention. 

In figure 10, a system of linear equations (represented by the matrix equation 
A o x = f Q ) y which have been derived from partial differential equations (PDEs) via 
finite difference quotients, can be characterized, for each iteration, as a seven step 
process. In figure 10, block 88, an initial matrix A Q and an initial residual vector 
^ 0 are obtained. Initial residual vector p 0 is the vector of current errors (residuals) 
for the system of nonlinear equations. Matrix A 0 is a sparse Jacobian matrix that 
describes the sensitivity of the residuals to changes in x for the system of nonlinear 
equations. In block 90 of figure 10, matrix y\ 0 is decomposed into smaller and 
smaller matrices until a matrix is found which is efficiently factored exactly. A 
matrix can be efficiently factored exactly when it has a bandwidth of less than or 
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equal to approximately 7. Starting with the largest decomposed matrix of A>that 
can be efficiently factored exactly (found in block 90 of figure 10), an approximate 
factorization for matrix ^ 0 is found such that A> = Mo = LU where L is a block 
lower triangular and U is a block upper triangular matrix, block 92 of figure 10. 

5 Using the LU factorization of matrix Ao (i.e., matrix Mo >' the value of A»* is 
computed where J is a computed approximation for the original linear system at the 
current time step, block 94 of figure 10. A predicted value of A>* is tnen 
computed for the next time step, block 96 of figure 10, which is then checked for 
convergence. If the predicted value Ao * converges (the'Ves'prong of block 98 in 

10 figure 10), the current time step has been completed, block 100 of figure 10. If the 
predicted value of Ao * does not converge (the"nd'prong of block 98), blocks 94 and 
96 of figure 10 are repeated. Blocks 88 through 100 of figure 10 are repeated for 
each time step in the simulation. 

15 In figure 1 1, a modified nested factorization technique to decompose matrix A>(U" 
block 92 of figure 10) is illustrated in figure 11. After initializing loop variable T to 
zero, block 102 of figure 11, a test is made to determine if the matrix denoted by A 
can be efficiently factored exactly. If matrix A cannot be efficiently factored 
exactly (the W prong of block 104 of figure 11), matrix A is partitioned (i.e. 
20 decomposed) such that A = L, + A +l + f/. where L; is a lower triangular 

matrix, JJM an upper triangular matrix, and A + i is a block diagonal matrix, block 
106 of figure 1 1 . Next, loop control variable T is incremented (block 1 08 in figure 
11) and processing continues at block 104 of figure 11. If matrix A can be 
efficiently factored exactly (the 'yes 1 prong of block 104 of figure 11), control 
25 variable V is assigned the current value of loop control variable T (block 1 10 in 
figure 1 1). For computational efficiency, matrix A. is preferably the largest (i.e., 
first found) decomposition of matrix Ao that can be efficiently factored exactly. 
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In figure 12, a technique for generating an approximate factorization of matrix ^4 0 
(block 92 in figure 10) is shown in figure 12. Matrix M „ is initialized to A»~& 
where matrix D is a diagonal matrix defined by: 

OP represents either a column-sum operator or a row-sum operator (block 1 12 of 
figure 12). The column-sum of matrix X is a diagonal matrix formed by summing 
the elements of matrix X in columns. The row-sum of matrix X is a diagonal matrix 
10 formed by summing the elements of matrix X in rows. Next, loop variable T is 
initialized to 'n-1' (block 1 16 of figure 12) and a test is made to determine if the 
factorization loop is complete. If variable T does not equal zero (the 'no 1 prong of 
block 1 18 of figure 12), the next M . matrix is computed in accordance with the 
following: 

15 

where I represent the identity matrix (block 120 of figure 12). Loop control variable 
T is then decremented (block 122), and processing continues at block 118. If 
20 variable T equals zero (the 'yes 1 prong of block 1 1 8), matrix M o < the factored 
approximation to A 0 ) is computed in accordance with: 

Mo=(LoM; , + ^(Af,+t/o). block 124offigure 12. 

25 In figure 13, during each iteration (step 94 through 100 of figure 10), an 

approximate solution x' to the original system of linear equations (A 0 * = T 0 ) is 
generated in accordance with figure 13. First, an initial solution s 0 is defined, e.g., 
S 0 = 0 (block 126 of figure 13). Next, loop variable T is initialized to 1 (block 128 
of figure 13) and a test is made to determine if the current solution has converged. If 
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the current solution (e.g., that solution associated with the current time step) has not 
converged (the 'no' prong of block 130), the equation MoSi = 7 m is solved for 
S t (block 132 of figure 13). One benefit of the methods illustrated in figures 1 1 and 
12 is that the following interim results are obtained at no additional computational 
5 cost: 

y = (LoMr 1 +/J ? 

Q = (Mi+Uq) '^and 
Zt-MtUtQ (fork=l„n-l) 

10 

These interim results can be used to compute a value for A J, » that is, 
Ao C M o 1 7 i-i) ' in accordance with the following (block 1 34 of figure 1 3): 

A(Mo 1 r.-,)=5 i +^ + Lo^-ZL t z t 

15 

The method of computing A J. outlined above and in figure 13 uses only 
(n + 1) diagonal vector multiplies, compared with the (2n + m) diagonal vector 
multiplies required by prior matrix multiplication techniques [where 'rf equals the 
value determined at block 1 10 of figure 1 1 , and W equals the number of diagonals 
20 (e.g., the bandwidth) of A n 1- The computational savings provided by the inventive 
method over the prior technique in two, three, and four dimensional (i.e., 2-D, 3-D, 
and 4-D respectively) oil and gas reservoir simulations can be described in 
percentage terms as shown in Table 1 below: 
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Table 1: Computational Savings of Inventive Concept 

m n %-Savings 



2-D Simulation 



40 



3-D Simulation -a 2 42.9 



4-D Simulation 



%-Savings = a:100 

2/i + m 



AAA 



In figure 13, block 136, a value for vector f , ™ a y be calculated using any suitable 
acceleration method. Illustrative stationary acceleration methods include the 
Jacobi, Guass-Seidel, successive overrelaxation (SOR), and symmetric successive 
10 overrelaxation (SSOR) methods. Illustrative nonstationary acceleration methods 
include the generalized minimal residual (GMRES), biconjugate gradient (BiCG), 
biconjugate gradient stabilized (Bi-CGSTAB), quasi-minimal residual (QMR), 
conjugate gradient squared (CGS), Chebyshev, and conjugate gradient (CG) 
techniques. 



If the value of f < generated in block 1 36 of figure 1 3 converges (the 'yes' prong of 
block 130 is active), the approximate solution x' associated with the current time 
step is computed in block 138 in figure 13 (i.e., Ao*' = 7 , is solved 
for x' ). Convergence may be defined in any suitable manner and, in one 
20 embodiment, is said to occur when the ratio between the magnitude of fi 

and r i+ i is less tnan 01% - Mathematically, this may be expressed in the following 



manner: 
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S4- < 0.001 , where fall represents the Euclidean norm of vector a . 
¥1 

Method steps of figures 10 through 13 may be performed by a computer processor 
executing instructions organized into program modules. Storage devices suitable for 
5 tangibly embodying computer program instructions include all forms of non-volatile 
memory including, but not limited to: magnetic disks (fixed and floppy); other 
magnetic media such as tape; optical media such as CD-ROM disks; and 

semiconductor memory devices such as EPROM and flash devices. 

10 The methods of this invention provide significant advantages over the current art. 
The invention has been described in connection with its preferred embodiments. 
However, it is not limited thereto. Changes, variations and modifications to the 
basic design may be made without departing from the inventive concepts in this 
invention. In addition, these changes, variations and modifications would be 
15 obvious to those skilled in the art having the benefit of the foregoing teachings. All 
such changes, variations and modifications are intended to be within the scope of 
this invention, which is limited only by the following claims. 
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WE CLAIM: 

L A simulation apparatus responsive to input data adapted for solving a 
system of nonlinear equations that represent a particular entity, said 
5 simulation apparatus generating a set of simulation results when said 
system of nonlinear equations are solved, said set of simulation results 
including one or more parameters which characterize said particular 
entity, comprising: 

10 first means for converting said system of nonlinear equations into a 
system of linear equations; 

second means for solving said system of linear equations thereby 
producing a set of results; and 

15 

third means for solving said system of nonlinear equations by using said 
set of results produced by said second means. 

2. The simulation apparatus of claim 1, wherein said input data includes 
20 reservoir data and said particular entity includes an oil or gas reservoir. 

3. The simulation apparatus of claim 1, wherein said second means for 
solving said system of linear equations comprises: 

25 means for obtaining a matrix Ao from the s y stem of linear equations; 

means for decomposing the matrix Ao until a factorable malrix Mo is 
obtained; 
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first computing means for computing an approximate factorization of the 
matrix Mo to produce a particular system of linear equations, which are 
represented by said matrix Mo- that approximates said matrix A,; and 

second computing means for computing a solution to said particular 
system of linear equations using said matrix Mo- 



4. The simulation apparatus of claim 3, wherein said third means solves 
said system of nonlinear equations by using said solution, computed by 

10 said second computing means, to said particular system of linear 
equations which use said matrix Mo- 

5. In a workstation apparatus adapted to receive reservoir data including 
a simulator adapted for receiving said reservoir data and for generating a 

15 set of simulation results which include one or more parameters 
representative of an oil or gas reservoir and a recorder or display 
apparatus adapted for recording or displaying said set of simulation 
results, said simulator comprising: 

20 first means for converting said system of nonlinear equations into 
a system of linear equations: 

second means for solving said system of linear equations thereby 
producing a set of results; and 

25 

third means for solving said system of nonlinear equations by 
using said set of results produced by said second means. 

6. The simulator of claim 5, wherein said second means for solving said 
30 system of linear equations comprises: 
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means for obtaining a matrix j\ Q from the system of linear equations; 

means for decomposing the matrix j\ Q until a factorable matrix M 0 is 
obtained; 

5 

first computing means for computing an approximate factorization of the 
matrix M 0 to produce a particular system of linear equations, which are 
represented by said matrix A/o » that approximates said matrix A 0 ; and 

10 second computing means for computing a solution to said particular 
system of linear equations using said matrix M o- 

7. The simulator of claim 6, wherein said third means solves said system 
of nonlinear equations by using said solution, computed by said second 

15 computing means, to said particular system of linear equations which use 
said matrix Mo- 

8. In a simulator, a method for solving a system of linear equations 
representing reservoir conditions, comprising the steps of: 

20 

(1) receiving a first signal representing physical characteristics of a 
reservoir; 

(2) obtaining, from the received signal, a residual vector f 0 representing 
25 errors associated with a system of nonlinear equations describing the 

reservoir; 

(3) obtaining, from the received signal, a first matrix A 0 representing 
sensitivity of the residual vector f : 0 to changes in the system of 

30 nonlinear equations describing the reservoir; 
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(4) decomposing matrix Ao in accordance with the following steps: 

(a) letting i = 0, 

(b) if matrix A can be efficiently factored exactly, letting n=i, 

5 else, 

(i) letting A = L, + A>, + £/; where L. « a lower block 
diagonal matrix and f/i is an upper block diagonal matrix, 

(ii) incrementing i, and 

(iii) repeating step (4)(b) for matrix A+i > 

10 (5) generating a second matrix Mo^at ^ an approximation to matrix 

Ao in accordance with the following steps: 

(a) letting M „ = A, ~ D where D is a dia S onal matrix defined 
by a specified operation on the matrix given by 

15 (b) generating, for i=(n-l) to 1, 

Mi = (Li+M »,>(' + Af £l/ «> • and 
(c) letting Mo=(LoA/r 1 + / )(M, + f/o)-' 

(6) generating, for k=l until convergence, the following values: 

20 (a) y = (LoM\ l +I) *r M . 

(b) q = (Mi+Uofy> 

(c) Zi=MnUi« (fori=l,...,n-l), 

(d) AoMon-^y + ^ + Lo^-S^z.-and 

(e) r M by accelerating the value A M o F t-i • 

25 
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(7) generating an approximate solution to the system of linear equations 



representing reservoir conditions using a value of said p t+1 that converged; and 

(8) generating a second signal representing physical characteristics of the 
5 reservoir. 

9. The method of claim 8, wherein the physical characteristics comprise 
reservoir pressure at one or more locations, water fraction at one or more 
locations, oil fraction at one or more locations, and gas fraction at one or 

10 more locations. 

10. The method of claim 8, wherein matrix A can be efficiently factored 
exactly when it has a bandwidth of less than or equal to approximately 7. 

15 11. The method of claim 8, wherein the specified operation comprises a 
column sum operation in accordance with the following: 



20 12. The method of claim 8, wherein the specified operation comprises a 
row sum operation in accordance with the following: 





25 13. The method of claim 8, wherein said accelerating comprises a 
generalized minimal residual method. 



14. A method for solving a system of linear equations representing 
reservoir conditions, comprising the steps of: 
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receiving a first signal representing physical characteristics of a reserv 

obtaining, from the received signal, a residual vector f 0 representing 
5 errors associated with a system of nonlinear equations describing the 
reservoir; 

obtaining, from the received signal, a first matrix A 0 representing 
sensitivity of the residual vector f o to changes in the system of 
10 nonlinear equations describing the reservoir; and 

decomposing matrix Ao in accordance with the following steps: 



(a) letting i = 0, 

(b) if matrix A can be efficiently factored exactly, letting n=i, 



15 



else, 




(i) letting A-L+Am + Ui where L, is a lower block 
diagonal matrix and f/^s an upper block diagonal matrix, 



(ii) incrementing i, and 

(iii) repeating step (4)(b) for matrix Am ■ 



15. The method of claim 14, further comprising the steps of: 



generating a second matrix Mo that is an approximation to matrix Ao 
in accordance with the following steps: 



25 



(a) letting M . - A." D wnere D is a dia S onal ™ atrix defined by a 
specified operation on the matrix given by 



2Xi-i M'iU 
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(b) generating, for i=(n-l) to 1, M. = (L, + M m^ 1 + MmV 
and (c) letting Mo = <Lo + /XA/.+W- 

16. The method of claim 15, further comprising the steps of: 
generating, for k=l until convergence, the following values: 

(a) y = (LoA/r 1 +/) Fw 

(b) * = (Mi+i/ (>)"?.• 

(c) Zi=MmUi« (fc>ri=l n-1), 

io (d) AoM^rw'y + ^ + Loff-lLzf.^ 

( e > F m b y accelerating the value A Q M o F *-i • 

17. The method of claim 14, wherein the physical characteristics 
comprise reservoir pressure at one or more locations, water fraction at 

15 one or more locations, oil fraction at one or more locations, and gas 
fraction at one or more locations. 

18. The method of claim 14, wherein matrix can be efficiently 
factored exactly when it has a bandwidth of less than or equal to 

20 approximately 7. 

19. The method of claim 14, wherein the specified operation comprises a 
column sum operation in accordance with the following: 

25 couuMytL^M'-U^ 
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20. The method of claim 14, wherein the specified operation comprises a 
row sum operation in accordance with the following: 



ROWSUM^ L,-, Ml t/i-,] • 



21. The method of claim 16 wherein said accelerating comprises a 
generalized minimal residual method. 

10 22. A computer program, residing on a computer readable medium, 
comprising instructions for causing a computer to: 



receive 



a first signal representing physical characteristics of a reservoir; 



15 obtain, from the received signal, a residual vector f „ representing errors 
associated with a system of nonlinear equations describing the reservoir; 

obtain, from the received signal, a first matrix Ac representing 
sensitivity of the residual vector r 0 to changes in the system of 
20 nonlinear equations describing the reservoir; and 

decompose said matrix A 0 » accordance with the following steps: 

(a) leti = 0, 

(b) if matrix Ai can be efficiently factored exactly, let n=i, 
25 else, 

(i) let A, = Li + A+i + U , wnere Li is a lower block 
diagonal matrix and J7, is an upper block diagonal matrix, 

(ii) increment i, and 

(iii) repeat step (4)(b) for matrix Am • 
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23. The computer program of claim 22, further comprising instructions 
to: 

5 generate a second matrix M 0 that is an approximation to matrix Ao * n 
accordance with the following steps: 

(a) letting M n = An~ D wh ere D is a diagonal matrix defined by 

n-l 

a specified operation on the matrix given by £ M^Um* 

(b) generating, for i=(n- 1 ) to 1 , A/, = (L, + M .♦, M+.M n U .) » 
10 and 

(c) letting Mo = (£oMr + / XM, + t/o)- 

24. The computer program of claim 23, further comprising instructions 
to: 



15 



generate, for k=l until convergence, the following values: 



(a) y = (L 0 Mr 1 +/)"r t -.. 

(b) ? = (Mi+£/ 0 )" , 5 i . 

20 (c) Zi=M?«UiV (fori a 1 n-1), 

(d) A Mo 1 r = y + ^ + L 0 ^ - E L. z. . and 

(e) r *+, by accelerating the value A 0 Mo P *-i • 
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